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1  Introduction 


Plasmas  have  the  potential  to  be  used  for  the  control  and  manipulation  of  diverse 
flow  related  applications.  Most  current  works  have  considered  the  prospect  of 
using  plasmas  for  flow  control  by  directly  affecting  the  mean  flow.  This  project 
is  studying  the  ability  of  plasmas  to  directly  influence  the  turbulent  fluctuations 
in  the  interior  of  flows.  Laser  energy  deposition  has  been  considered  to  be  the 
mode  of  plasma  generation. 

Laser  energy  deposition  has  been  studied  in  great  detail  (eg.  Maker  et  al 
1963,  Damon  k  Tomlinson  1963,  Meyerand  k  Haught  1963).  Applications  of 
laser-induced  breakdown  of  a  gas  include  localized  flow  control  of  supersonic 
flows  ( Adelgren  et  al  2003) ,  drag  reduction  in  supersonic  and  hypersonic  flows 
(Riggins  1999),  ignition  of  combustion  gases  (Phuoc  2000)  and  providing  thrust 
to  aerospace  vehicles  (Molina-Morales  et  al  2001,  Wang  et  al  2001). 

When  a  laser  beam  is  focused  on  a  small  volume  of  gas,  the  gas  molecules  in 
the  focal  volume  absorb  energy  and  get  ionized.  The  laser  beam  can  convert  an 
optically  transparent,  gaseous  region  into  opaque,  highly  conducting  plasma,  in 
a  few  nano-seconds.  The  plasma  formation  process  involves  multiple  steps.  The 
initial  release  of  electrons  is  due  to  multi-photon  ionization,  which  requires  that 
multiple  photons  are  incident  on  an  atom  simultaneously.  Multiple  photons  are 
needed  since  ionization  energies  for  most  gases  are  larger  than  the  energy  of 
a  single  photon  in  the  laser  beam.  During  this  process,  the  electron  number 
density  increases  linearly  in  time.  The  released  electrons  in  the  focal  volume 
absorb  laser  energy  due  to  inverse  bremsstrahlung  absorption,  where  a  free 
electron  in  presence  of  a  third  body,  absorbs  energy  and  gets  excited.  Many  such 
interactions  are  needed  for  the  electron  to  gain  sufficient  energy  to  impact  ionize 
neutral  atoms.  The  resultant  electron  concentration  increases  exponentially  in 
time. 

The  plasma  thus  formed,  absorbs  laser  energy  and  spreads  out  in  the  direction 
of  the  laser  beam  due  to  a  combination  of  absorption  and  reflection.  The  highly 
ionized  plasma  in  the  focal  volume  reflects  part  of  the  laser  energy  incident  on  it. 
This  energy  is  absorbed  by  adjacent  molecules,  in  the  direction  of  the  laser  beam. 
These  molecules  get  ionized,  and  start  reflecting  laser  radiation.  This  process 
continues  until  the  plasma  has  spread  out  and  assumed  a  tear-drop  shape.  The 
energized  electrons  collide  with  the  heavy  particles  leading  to  heating  of  the 
gas.  Also  due  to  these  collisions,  the  electrons  get  de-energized  and  recombine 
with  ions,  sharply  decreasing  the  electron  number  density.  Thus  an  energy  spot 
with  temperature  and  pressure  higher  than  that  of  the  surroundings  is  left  at 
the  end  of  the  plasma  formation  process.  The  sharp  pressure  gradients  lead  to 
formation  of  a  blast  wave  that  propagates  into  the  background  gas. 

Recent  experimental  work  on  pulsed  laser  induced  breakdown  in  quiescent  air  or 
nitrogen,  include  Jiang  et  al.  (1998),  Lewis  et  al.  (1999),  Adelgren  et  al.  (2001) 
and  Yan  et  al.  (2003).  Jiang  et  al  (1998)  used  a  laser  beam  of  L38  J  focused 
on  a  3  mm  diameter  spherical  region  to  cause  breakdown  in  air.  The  laser  was 
pulsed  for  a  duration  of  18  nanoseconds.  Adelgren  et  al.  (2001)  used  a  Nd;YAG 
laser  of  200  mJ  pulsed  for  10  nanosecond  duration  in  air.  Yan  et  al.  (2003) 
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LASER 


Figure  1:  Schematic  illustrating  the  thermal  interaction  between  isotropic  tur¬ 
bulence  and  spherical  and  conical  idealizations  of  a  region  of  plasma. 


used  a  108  -  145  mJ  Nd;YAG  laser  pulsed  for  a  duration  of  10  nanoseconds 
over  a  focal  volume  of  3  mm3  to  induce  breakdown  in  air.  The  experimental 
data  show  a  wide  separation  in  time-scales  between  the  laser  pulse  duration 
and  observations  of  the  blast  wave  propagation.  As  seen  above ,  the  laser  is 
pulsed  on  a  time-scale  of  10  nanoseconds.  The  blast  wave,  on  the  other  hand  is 
observed  on  a  time-scale  of  10  to  100  microseconds.  Since  plasma  formation  is 
on  a  time-scale  of  the  same  order  as  the  laser  pulse  duration,  the  three  to  four 
order  of  magnitude  separation  in  time-scale  suggests  that  the  plasma  may  be 
considered  to  be  formed  instantaneously  to  a  first  approximation,  to  evaluate 
its  gas-dynamic  effect  on  the  surrounding  fluid. 

Various  theoretical  models  have  been  used  to  understand  different  features  of 
this  phenomenon.  Similarity  laws  proposed  by  Taylor  (1950),  Sedov  (1959)  and 
Von  Neumann  (1963)  for  the  radius  and  speed  of  the  blast  wave  produced  from 
a  point  source  apply  to  the  early  stages  of  the  blast  wave,  when  it  is  strong. 
Note  however  that  at  early  stages,  the  blast  wave  need  not  be  spherical,  which 
reduces  applicability  of  the  classical  analysis.  Erode  (1955)  performed  a  numer¬ 
ical  simulation  of  the  blast  wave  and  concluded  that  the  ideal  gas  assumption 
was  reasonable  for  shock  pressures  less  than  10  atmospheres  in  air.  Steiner  et 
al,  (1998)  perform  computations  using  a  real  gas  model  to  show  that,  when 
initialized  with  a  self  similar  strong  shock  solution,  the  shock  radius  in  the  real 
gas  model  is  quite  close  to  that  predicted  by  the  classical  point  source  explosion 
in  an  ideal  gas.  Other  computations  of  blast  wave  propagation  in  quiescent  air 
include  those  by  Jiang  et  al  (1998)  and  Yan  et  al.  (2003).  Dors  et  al.  (2000) 
present  a  computational  model  which  considers  the  asymmetry  of  laser  energy 
deposition  as  well  as  ionization  and  dissociation  effects  on  fluid  properties.  The 
initial  stages  of  plasma  formation  due  to  laser  energy  deposition  were  modeled 
by  Kandala  and  Candler  (2002). 

This  work  considers  generation  of  laser  induced  plasma  in  quiescent  air  and  in 
background  isotropic  turbulence,  A  schematic  for  the  problem  is  shown  in  figure 
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1 .  Local  thermodynamic  equilibrium  conditions  are  assumed  to  apply.  Only  the 
thermal  effects  of  the  plasma  are  considered.  Two  types  of  idealizations  have 
been  used  in  the  simulations  to  represent  the  initial  shape  of  the  plasma  - 
spherical  and  tear-drop  shaped.  The  spherical  idealization  corresponds  to  a 
point  plasma  region  and  is  a  good  approximation  only  if  the  size  of  the  plasma 
is  small  compared  to  other  length-scales  in  the  flow.  The  tear-drop  shaped 
idealization  on  the  other  hand  resembles  the  initial  shape  of  the  plasma  observed 
in  experiments  of  laser  induced  break  down  of  a  gas, 

A  Fourier  spectral  solver  was  developed  for  the  simulations.  The  solver  assumes 
periodic  boundary  conditions.  If  the  plasma  is  assumed  to  be  spherical  in  nature, 
then  the  flow  would  vary  radially  and  in  time.  If  the  plasma  is  assumed  to  be 
tear-drop  shaped,  the  flow  would  be  axisymmetric  and  non -station ary  in  time. 
For  both  cases,  the  periodic  boundary  conditions  would  work  fine  as  long  as  the 
blast  wave  does  not  touch  the  domain  boundaries.  This  is  because  if  the  blast 
wave  is  well  resolved,  the  gradients  ahead  of  it  will  be  zero. 

This  grant  was  used  to  support  the  doctoral  research  of  Mr.  Shankar  Ghosh, 
Publications  associated  with  this  work  are  listed  below. 

•  Direct  numerical  simulation  of  the  thermal  effects  of  plasmas  on  turbulent 
flows  {AIAA  Paper  2005-0407), 

•  Direct  Numerical  Simulation  of  the  thermal  effects  in  plasma /turbulence 
interaction  (AIAA  Paper  2007-993). 

•  DNS  of  the  thermal  effects  of  laser  energy  deposition  in  a  fluid  at  rest 
(submitted  to  Journal  of  Fluid  Mechanics}. 

•  Numerical  issues  associated  with  simulation  of  strong  blast  waves  in  high 
temperature  flows  interacting  with  a  turbulent  background  (under  prepa¬ 
ration  for  submission  to  the  Journal  of  Computational  Physics}. 

The  principal  contributions  of  this  work  are  listed  below. 

•  A  parallel  compressible  Navier-Stokes  solver  has  been  developed.  The 
solver  uses  Fourier  spectral  spatial  derivatives.  SkewT  symmetric  form  of 
the  non-linear  terms  are  used  to  suppress  aliasing  error,  A  fourth  order 
Runge-Kutta  scheme  has  been  used  for  time  advancement. 

•  A  pencil  data  structure  has  been  used  to  parallelize  the  solver.  An  efficient 
algorithm  has  been  implemented  to  compute  transposes  for  data  along  the 
y  and  z  directions.  The  solver  has  been  ported  and  tested  on  different 
platforms  (Sdsc  Teragrid,  Datastar,  Bluegene),  the  largest  simulation  so 
far  being  a  512  cube  simulation  on  1024  processors, 

•  A  p red ictor-cor rector  based  shock  capturing  scheme  has  been  implemented 
to  capture  strong  shock  waves  observed  in  the  simulations.  The  scheme 
(Yee  et  al  1999}  has  a  finite  volume  based  formulation  and  has  been  applied 
as  a  corrector  step  to  the  spectral  base  numerical  scheme. 
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•  Shock  capturing  needs  to  be  implemented  without  excessively  dissipating 
the  turbulence  in  the  background.  This  is  achieved  by  developing  a  non¬ 
linear  switch  that  detects  the  shock  wave  and  applies  shock  capturing  only 
to  its  vicinity, 

•  Since  the  simulations  involve  very  high  temperatures,  the  finite  volume 
formulation  of  the  shock  capturing  scheme  needs  to  be  extended  suitably 
for  high  temperatures.  A  formal  derivation  of  the  Jacobian  matrix,  eigen 
values  and  eigen  vector  matrices  has  been  derived  for  a  generalized  geom¬ 
etry  (Appendix). 

•  There  is  significant  expansion  in  the  plasma  core  leading  to  very  low  densb 
ties  there.  Slight  oscillations  can  thus  lead  to  densities  becoming  negative 
rendering  the  solution  unstable.  To  avoid  this  a  logarithm  formulation 
of  the  continuity  equation  has  been  used  to  ensure  that  density  remains 
positive  even  in  regions  of  very  low  density. 

•  Based  on  the  amount  of  energy  deposited  in  the  flow,  the  blast  wave 
reaches  the  domain  boundaries  on  a  certain  time  scale.  To  study  further 
evolution  of  the  plasma  core  the  domain  boundaries  need  to  be  expanded. 
This  is  done  by  interpolating  the  solution  form  a  smaller  to  a  larger  domain 
using  a  spectral  interpolation  scheme. 

•  In  the  case  of  spherical  energy  deposition,  the  post  energy  deposition  pro¬ 
cess  can  be  divided  into  two  parts  -  shock  formation  and  shock  propaga¬ 
tion.  These  phases  have  been  described  based  on  conversion  of  internal 
energy  of  the  core  into  kinetic  energy  of  surrounding  fluid  elements  and 
generation  of  sharp  pressure  gradients  in  the  flow. 

•  The  shock  radius  evolution  in  time  has  been  compared  to  the  strong  shock 
solution  (Taylor  1950,  Sedov  1959)  and  found  to  compare  well  at  short 
times  and  deviate  towards  the  acoustic  limit  at  longer  times.  Lesser  the 
energy  deposited  in  the  flow,  sooner  is  the  deviation  towards  the  acoustic 
limit.  Similar  analysis  is  also  presented  for  velocity  at  the  shock  front. 

•  A  reverse  flow  is  observed  behind  the  shock  front  at  later  times.  Existence 
of  this  reverse  flow  is  explained  based  on  conservation  of  mass  behind  the 
shock  front.  The  reverse  flow  is  observed  to  be  symmetric  in  nature, 

•  In  the  case  of  the  tear-drop  shaped  plasma,  shock  radius  evolution  in 
time  and  jumps  in  flow  variables  at  the  shock  front  are  compared  with  data 
obtained  from  experiments  conducted  by  Adelgren  (2001)  and  simulations 
conducted  by  Kandala  and  Candler  (2003). 

■  Again  post  energy  deposition  is  divided  into  formation  and  propagation  of 
a  blast  wave.  However  the  blast  wave  is  initially  asymmetrical  tear-drop 
shaped  as  characterized  by  sharp  pressure  gradients  normal  to  the  plasma 
axis.  This  blast  wave  propagates  into  the  background  becoming  spherical 
in  time. 
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•  Barodinic  vorticity  is  found  to  be  generated  near  the  leading  and  trailing 
edges  of  the  plasma  at  short  times-  This  is  explained  based  on  a  relation 
suggested  by  Truesdeli  (1952). 

•  Budgets  have  been  computed  for  the  vorticity  transport  equation.  They 
indicate  that  vorticity  generated  by  barodinic  means  is  convected  radially 
outward  and  suppressed  by  bulk  expansion  behind  the  blast  wave. 

•  For  longer  times  an  asymmetric  reverse  flow  is  observed  behind  the  shock 
front  resulting  in  the  rolling  up  of  the  flow  field.  This  leads  to  advec- 
tion  of  the  temperature  profile  such  that  the  maximum  temperature  shifts 
away  from  the  plasma  axis.  This  is  in  good  agreement  with  experimental 
observation  of  the  resulting  flow  field, 

•  The  resulting  vorticity  field  is  observed  to  be  much  more  complicated  than 
a  simple  toroidal  vortex  ring  observed  in  the  experiments.  However  the 
vorticity  magnitude  is  observed  to  be  much  larger  for  the  rings  closest 
to  the  plasma  axis  and  so  these  are  the  ones  that  are  observed  in  the 
experiments. 

•  In  the  case  of  the  spherical  plasma  interacting  with  background  turbu¬ 
lence,  a  spherically  symmetric  blast  wave  is  observed  to  propagate  through 
the  background  turbulence.  Turbulence  levels  are  observed  to  get  ampli¬ 
fied  across  the  blast  wave.  Amplification  is  most  pronounced  in  the  radial 
component  of  velocity  and  the  transverse  component  of  vorticity. 

•  In  the  case  of  the  tear-drop  shaped  plasma  interacting  with  background 
turbulence,  a  tear-drop  shaped  blast  wave  is  observed  to  propagate  through 
the  background  turbulence  becoming  more  and  more  spherical  in  time. 

•  The  blast  wave  strongly  compresses  the  background  turbulence.  Signifi¬ 
cant  expansion  is  observed  in  the  turbulent  core. 

•  The  blast  wave  is  initially  very  strong  and  is  not  distorted  by  the  back¬ 
ground  turbulence.  However,  as  it  propagates  into  the  background,  the 
effect  of  the  turbulence  is  to  dampen  its  intensity  and  soon  it  becomes 
weak  enough  to  get  distorted  through  interaction  with  the  turbulence. 

•  Statistics  computed  for  turbulent  kinetic  energy'  and  divergence  of  velocity 
fluctuations  indicate  that  turbulence  levels  get  amplified  in  the  vicinity  of 
the  blast  wave  due  to  local  compression  and  get  suppressed  in  the  core 
due  to  expansion  there.  This  trend  is  supported  by  linear  analysis  of 
turbulence  going  through  a  shock  wave  (Mahesh  et  al.  1995,1997). 

•  Statistics  for  vorticity  fluctuations  indicate  amplification  in  vorticity  levels 
near  the  leading  edge  of  the  plasma  where  the  mean  vorticity  levels  were 
observed  to  be  high. 
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This  report  is  organized  as  follows.  Section  2  describes  details  of  the  algorithm 
used,  validation  tests  performed  and  numerical  issues  addressed.  Section  3  con¬ 
tains  simulation  results*  Section  3*1  describes  the  spherical  idealization  of  the 
plasma  while  section  3.2  contains  description  of  the  tear-drop  shaped  idealiza¬ 
tion  of  the  plasma*  Section  3.3  describes  interaction  of  the  spherical  idealization 
of  the  plasma  with  background  isotropic  turbulence*  Section  3*4  contains  de¬ 
scription  of  the  tear-drop  shaped  idealization  of  the  plasma  interacting  with 
isotropic  turbulence  in  the  background*  A  short  summary  is  presented  in  sec¬ 
tion  4. 


2  Simulation  methodology 

2.1  Algorithm  description 

The  governing  equations  are  the  continuity,  and  compressible  Navier-Stokes 
equations  applied  to  a  real  gas, 


dp  dpuj  _ 

at  +  dxi 


1 

dpUjUi 

d_ 

pSij  — 

P(T)  ( 

duj 

dt 

dxi 

dxi 

Re  \ 

^dxj 

dx% 

dpeT 

dpe  rtii 

=  —  \ 

-pUj  + 

p{T)  ( 

'fhH 

diij 

dt 

dxi 

dxi  [ 

Re  ( 

,dx} 

a 7, 

9  ( 

h(T) 

dT 

\ 

dxi  \  (7  -  l)RePr  / 


(1) 


2  duk 
3dTk°ij 
2  du^  r 

3dxk  ,J 


(2) 

-K3) 


where  all  the  variables  are  non-dimensional  ized  with  respect  to  their  initial 
background  values. 


Xi  =  Xi‘/Lo’,  Uj=ttj*/co*,  t  =  t‘co‘/Lo',  (4) 

P  =  p7po*.  p  =  p7po*co*2,  r  =  T*/T0* 
n(T)  =  p(r)7po*,  =  «(r)7«0'. 

Here,  the  subscript  *0’  denotes  initial  background  values  and  the  superscript,  *** 
denotes  dimensional  variables.  Lq  is  the  length  scale  and  is  obtained  by  com¬ 
paring  the  non-dimensional  length  of  the  plasma  region  used  in  the  simulations 
to  the  actual  dimensional  length  of  the  plasma,  cj  is  the  speed  of  sound  based 
on  initial  background  temperature;  i.e. 

CO*  =  (70  W),/2.  (5) 

and  k(T)¥  are  dimensional  coefficients  of  viscosity  and  thermal  conduc¬ 
tivity  obtained  by  assuming  an  equilibrium  model  for  air  (figure  2), 

The  equation  of  state  is 
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Figure  2;  (a)  Variation  of  the  coefficient  of  viscosity  \i  with  temperature  7\  (b) 
Variation  of  the  coefficient  of  thermal  conductivity  k  with  temperature  I\ 


p  =  PR{T )T,  {6) 

where 

R(T)  =  R'(T)h ofio,  (7) 

and  the  variation  of  R*  with  temperature  is  shown  in  figure  3b.  The  total  energy 
is  related  to  internal  energy  and  kinetic  energy  as 

per  -  pe^-^pUiUi.  (8) 

and  temperature  is  obtained  form  internal  energy  using  the  equilibrium  depen¬ 
dence  of  internal  energy  on  temperature  shown  in  figure  3a. 

The  Reynolds  number  and  Prandtl  number  are  given  by 


Re^  po* cq'Lq* /pQ\  Ft  ^  ^0*Cpo7«o‘*  (9) 

Cpo*  is  the  specific  heat  at  constant  pressure  at  T *  =  TqV 

The  Navier-Stokes  equations  are  solved  using  Fourier  methods  to  compute  the 
spatial  derivatives.  Any  variable  /  is  discretely  represented  as 


N  i/2-1  Na/ 2-1  N3/2-i 

/(*  1,12,13)*  Y  Y  Y  f(kuk3,k3)e<k'x'+k’x’+k>xJ  (10) 

*,  =  — Ni/2  Ar2=-N3/2Jfc*^-/V>/2 
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e’/Co 


(a)  (b) 


T*/T0m  T*/Tq* 

Figure  3:  (a)  Variation  of  internal  energy  e  with  temperature  7\  (b)  Variation 
of  R  with  temperature  T\ 


where  /( A;1t  ^2,^3)  are  the  Fourier  coefficients  of  /,  and  A/yt  N2  and  /V3  are  the 
number  of  points  used  to  discretize  the  domain  along  Xir  X2  and  X3  respectively. 
The  Fourier  coefficients  of  the  spatial  derivatives  are  therefore 


bL-uuT  JiZ- 

—  lKaj T 

(JXft  C/X^X  a 


(11) 


A  collocated  approach  is  used,  and  the  solution  is  advanced  in  time  using  a 
fourth  order  Runge-Kutta  scheme.  The  skew-symmetric  form  of  the  convection 
terms 


?Il+fJ*L+0*l 

dxj  J  dxj  ^  dxj 


(12) 


9/g  =  1 

dxj  2 

is  used  to  suppress  aliasing  errors  resulting  from  the  nonlinear  convection  terms 
(Blaisdel)  1991).  The  above  algorithm  is  implemented  for  parallel  platforms 
using  MPI.  The  library  FFTW  is  used  to  compute  Fourier  transforms,  and  a 
pencil  data  structure  is  used.  Each  processor  stores  data  along  the  entire  extent 
of  the  x\  —  direction,  while  data  along  the  X2—  and  £3—  directions  are  equally 
distributed  among  the  processors,  Fourier  transforms  along  the  Xj  -  direction 
are  therefore  readily  computed,  while  transforms  in  the  other  directions  require 
that  the  data  he  transposed  prior  to  transforming. 


2.2  Shock  capturing 

Recall  that  a  strong  shock  wave  propagates  through  the  flow  domain,  when 
energy  is  added  instantaneously.  Experiments  in  laser  induced  breakdown  (Yan 
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et  al  2003)  show  that  the  maximum  temperature  in  the  core  is  very  high.  This 
leads  to  sharp  gradients  in  the  flow  variables.  Since  the  flow  solver  uses  spectral 
methods  for  spatial  discretization,  resolving  these  sharp  gradients  requires  a 
highly  refined  mesh.  The  computational  cost  therefore  increases  significantly 
with  increasing  core  temperatures.  The  Fourier  spectral  method  is  therefore 
combined  with  a  shock  capturing  scheme  proposed  by  Yee  et  al  (1999),  to  avoid 
resolution  of  the  shock  thickness. 

The  shock  capturing  scheme  is  based  on  the  finite  volume  methodology,  and 
is  applied  as  a  corrector  step  to  the  Fourier  discretization  used  in  this  work. 
In  the  first  step,  the  predicted  form  of  the  solution  vector  is  obtained  using 
Fourier  methods  as  discussed  in  the  previous  section*  This  solution  vector  is 
then  corrected  using  the  filter  numerical  fluxes  obtained  by  using  a  characteristic 
based  filter 


I/n+l  =  +  Af 


(13) 


(^ijVA+1/2 


The  filter  numerical  flux  vector  is  of  the  form 

where  R  is  the  right  eigen  vector  matrix.  The  elements  of  4>*  are  denoted  by 
4>u  and  are  given  by 


“  *^t+l/2Tj,fc^!+l/2oJfc  (15) 

The  parameter  k  is  problem  dependent  and  lies  between  0,03  and  2  (Yee  et  al 
1999).  x  “  1.0  has  been  used  in  the  simulations.  The  function  jk  ls 
Harten  switch  (Harten  1978)  and  depends  on  the  Left  eigen  vector  matrix  L. 
The  formulation  used  for  6*ven  by  the  Harten- Yee  upwind  TVD 

form  (Yee  et  al  1999), 

However  for  high  temperature  flows,  separate  treatment  needs  to  be  provided 
to  the  computation  of  the  eigen  vector  matrices  R  and  L.  The  specific  heats  at 
constant  pressure  and  volume  Cp  and  Cv  are  no  longer  constants  but  depend 
strongly  on  temperature  (figure  4).  Also  the  internal  energy  and  enthalpy  can 
no  longer  be  obtained  from  the  simple  relations 

e  —  CyT  \  h  =  CPT  (16) 

The  sound  speed  is  no  longer  obtained  by 

c  =  (7flr)1/2  (17) 
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Cp'fy oRo 


(b) 


T'/To* 


Figure  4:  (a)  Variation  of  Cp  with  temperature  7\  (b)  Variation  of  Cy  with 
temperature  T. 


and  R  no  longer  remains  a  constant  but  becomes  a  function  of  temperature 
(figure  3b).  Thus  to  compute  the  eigen  vector  matrices,  the  Jacobian  matrix 
^  needs  to  be  recomputed.  Here  F  denotes  the  flux  vector  and  U  denotes  the 
vector  of  the  conserved  variables. 

Let  the  equation  of  state  for  the  gas  be  given  by 


p  =  pR(T)T .  (18) 

Then  the  elemental  change  in  pressure  dp  can  be  written  as 

dp=  RTdp  +  p^R  +  T^,)dT  (19) 

Now  since  e  =  e(T)  only 

Where  the  internal  energy  e  can  be  written  in  terms  of  the  conservative  variables 
Uj.  Using  the  above  relations  the  Jacobian  matrix  along  the  x-direction  can  be 
computed  as 


0 

1 

0 

0 

0 

\ 

(-«2  +  RT  -  A(e  -  ek)) 

(2  -  4)  *  ti 

~Av 

-Aw 

A 

Jx  = 

—vu 

V 

u 

0 

0 

—wu 

w 

0 

u 

0 

-(e0  +  A(e-  e*))u 

(eo+RT-Au2) 

-Auv 

—  Auw 

(1  +  A)u 

) 

Where 

e0«fi  +  cjt  (21) 
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is  the  total  energy  and  ek  is  the  kinetic  energy  given  by 


1 

ek  =  - UiUi 


(22) 


And  the  factor  A  is  given  by 


(23) 


Similar  Jacobian  matrices  can  be  computed  along  y  and  z  directions  as  well. 
Note  that  evaluation  of  these  Jacobian  matrices  need  knowledge  of  the  variation 
of  e  and  R  with  temperature.  This  can  be  obtained  for  air  assuming  it  to  be 
composed  of  several  species  present  under  conditions  of  chemical  and  thermal 
equilibrium  (figure  3). 

The  eigen  values  of  the  above  Jacobian  matrix  is  given  by 


A  =  (  u  —  Ci ,  u  u  +  Ci ,  it  u  ) 
Where  the  modified  speed  of  sound  cx  is  given  by 


c  =  ((i  +  A)i?r),/a 


(24) 


it  may  be  noted  that  C\  reduces  back  to  c  for  the  low  temperature  conditions 
and  the  Jacobian  matrix  reduces  back  to  the  low  temperature  Jacobian.  The 
modified  sound  speed  depends  on  the  variation  of  e  and  R  with  T  through 


the  factor  A ■  c3  has  been  computed  and  compared  to  data  for  speed  of  sound 


obtained  assuming  air  as  a  mixture  of  multiple  species  in  equilibrium  (figure  5). 
A  reasonable  comparison  is  obtained  up  to  T*  —  15000 A'. 

The  right  eigen  vector  matrix  R  can  be  obtained  by  solving  the  system  of 
equations 


[J]Ri  =  A,  A, 


(25) 


Where  A{  are  the  individual  eigen  values  and  A,  are  the  eigen  vectors.  Define  a 
set  of  variables  ej  and  £2  such  that 


(26) 


And 


€2  —  e  -  e\ 


(27) 


Also  define  e**  such  that 


cjt  =  e*  —  C2 

Then  the  right  eigen  vector  matrix  for  the  x  direction  is  obtained  as 


(28) 
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1 

1 

1 

0 

0  \ 

U  -  Cl 

u 

U  +  Ci 

0 

0 

V 

V 

V 

-1 

0 

w 

w 

w 

0 

I 

(ho  -  Ciu) 

ej t* 

(ho  +  cju) 

-v 

w  / 

Where  ho  is  the  total  enthalpy  given  by 

h0  =  h  +  e*. 

Then  the  left  eigen  vector  matrix  is  obtained  as 

\D  =  [i*r‘ 

and  for  the  x-direction  it  is  given  by 


/ 

_  Av 

Aw 

A 

f 

2cj  3 

~ zpr 

2ci 7 

2^ 

2c  i7 

Ci2  ^ 

■4ti 

Av 

Aw 

A 

*  CJ2 

fi7 

cl7 

cT7 

Ar.k  ~  ciu 

Av-ci 

_  Av 

_  Aw 

A 

2C|* 

2 ci7 

2c  i 2 

2c  i 2 

V 

0 

-I 

0 

0 

\ 

0 

a 

i 

0 

(29) 


(30) 


Similarly  eigen  vectors  Ry,  Rz ,  Ly  and  Lz  can  be  computed  along  the  y  and  z 
directions  respectively.  A  more  generalized  derivation  that  applies  to  generalized 
geometries  is  presented  in  the  Appendix. 


2,3  Logarithm  formulation  of  the  continuity  equation 

When  laser  energy  is  added  to  a  flow  at  rest,  there  is  noticeable  expansion 
at  the  core.  This  results  in  very  small  values  of  density  at  the  core.  When 
the  continuity  equation  was  advanced  in  time  with  density  as  the  dependent 
variable,  the  solution  was  Found  to  become  unstable.  It  was  therefore  decided 
to  solve  for  the  logarithm  of  density  as  the  variable.  Define 


u  =  Iogp  =>  p  =  ev  (31 ) 

The  continuity  equation  becomes 

dv  dv  _  dui 
dt  +  Utdx%  dx%  * 

Note  that  p  is  always  positive  when  computed  as  ev,  even  for  very  small  values  of 
p.  The  logp  formulation  of  the  continuity  equation  therefore  makes  the  solution 
stable  in  regions  of  very  low  density. 
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T* 


Figure  5:  (a)  Comparison  of  c\  with  data  for  speed  of  sound,  - :  speed 

of  sound  obtained  from  data,  ■  :  computed  speed  of  sound  — - :  c  = 

{■yoRoTf2 

2.4  Validation  Results 

The  solver  with  shock  capturing  and  the  logp  formulation  for  the  continuity 
equation  is  validated  using  a  periodic  shock  tube  problem  and  isotropic  turbu¬ 
lence.  Figure  6a  shows  instantaneous  profiles  of  the  solution  for  the  periodic 
shock  tube  problem.  A  region  of  high  density  fluid  is  located  at  the  center  of 
a  tube.  This  region  is  separated  from  low  density  regions  of  the  same  fluid  on 
either  side.  The  initial  density  jump  is  from  0.3  to  3  and  viscosity  in  the  problem 
is  0.005.  The  initial  pressure  profile  is  obtained  from  the  ideal  gas  law,  keeping 
temperature  constant.  This  initial  condition  generates  shock  waves,  which  prop¬ 
agate  in  either  direction.  Density,  pressure  and  temperature  profiles  obtained  at 
t  =  1.0  using  shock  capturing  and  the  log  p  formulation  are  compared  to  those 
obtained  without  using  these  formulations.  The  profiles  show  good  agreement. 
Figure  fib  shows  decay  of  turbulent  kinetic  energy  for  isotropic  turbulence.  The 
initial  velocity  fine tuat ions  are  generated  using  Rogallo’s  (1981)  method.  The 
initial  energy  spectrum  is  top- hat,  with  energy  in  wave  numbers  between  8  and 
16.  The  turbulent  Reynolds  number  is  50  and  the  fluctuation  Mach  number  is 
0.3.  Note  that  decay  of  turbulent  kinetic  energy  in  time  with,  and  without  the 
logp  formulation  of  the  continuity  equation  compare  well. 

2.5  Validity  of  Navier-Stokes  simulation 

This  work  considers  only  the  thermal  effects  of  the  plasma,  assuming  that  these 
effects  are  decoupled  from  the  effects  of  charged  particles  in  the  How,  This  can 
be  justified  by  recalling  that  the  time- sc  ale  for  plasma  generation  is  significantly 
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Figure  6:  {a)  Comparison  of  profiles  for  the  periodic  shock  tube  problem  at 

t  s=s  LG,  -  without  shock  capturing  and  logp  formulation  ■  with 

shock  capturing  and  logp  formulation,  (b)  Decay  of  turbulent  kinetic  energy 
for  isotropic  turbulence.  Time  is  non-dimensional  ized  by  the  eddy  turn  over 
time  r,  -  without  logy?  formulation,  ■  with  log/?  formulation. 


smaller  than  the  time-scale  on  which  the  gas-dynamic  effects  on  the  flow  are 
observed.  Also  at  the  end  of  plasma  formation,  the  electron  number  density 
decreases  sharply  due  to  electron  recombination.  Typically  for  ND-YAG  laser 
plasma  at  an  ambient  pressure  of  1  atmosphere,  the  electron  number  density 
decreases  to  less  than  one  tenth  of  it's  maximum  value  in  100  nanoseconds.  So  at 
the  end  of  plasma  formation,  all  that  is  left  is  a  region  of  high  temperature  and 
pressure.  If  the  simulations  are  initialized  with  the  profiles  of  thermodynamic 
quantities  at  the  end  of  the  plasma  formation  process,  then  the  effects  of  charged 
particles  will  not  be  very  significant. 

Note  that  the  temperature  in  the  core  of  the  plasma  is  quite  high  at  initial  times. 
This  leads  to  finite  rate  chemical  reactions  in  the  core  of  the  plasma.  Heat  ab¬ 
sorbed  or  released  during  these  reactions  affects  the  thermodynamic  properties 
in  the  core.  However,  the  core  cools  down  fast.  Also,  temperatures  near  the 
blast  wave  are  much  lower  compared  to  those  in  the  core  of  the  plasma.  Thus  for 
the  purpose  of  studying  the  properties  of  the  blast  wave,  the  effects  of  chemistry 
would  not  be  significant.  This  will  be  demonstrated  by  comparing  shock  radius 
and  jumps  in  the  thermodynamic  quantities  to  experimental  results.  Also,  the 
ideal  gas  law  has  been  assumed,  and  fluid  properties  like  viscosity  and  ther¬ 
mal  conductivity  have  been  assumed  to  vary  as  a  power  law  with  temperature. 
Again,  these  assumptions  hold  after  an  initial  time  when  the  plasma  generation 
is  complete.  Local  thermodynamic  equilibrium  conditions  have  been  assumed. 
These  conditions  are  reasonably  valid  after  a  few  nanoseconds  owing  to  the  high 
temperatures  involved. 
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Figure  7:  (a)  Contours  of  pressure  showing  propagation  of  blast  wave  in  flow  at 
rest,  (b)  Illustration  of  grid  convergence.  Radial  profiles  of  density  for  different 

mesh  spacing,  - :  Ax  =  0.0157,  ■  :  Ar  —  0.0196,  — --  ;  Ax  =  0.0245, 

- :  Ax  =  0.0491. 

3  Results 

3.1  Spherical  idealization 

A  spherical  region  of  the  flow  is  heated  at  the  center  of  the  domain.  The  initial 
temperature  distribution  has  a  Gaussian  profile: 

T  -_ 1  =  {T0  -  l)e-r,/r»  (33) 

AT  ATo 

where  the  temperatures  are  non -dimensional!  zed  by  the  ambient  temperature. 
AT  is  the  local  temperature  excess  over  the  ambient  value,  and  AT0  is  the 
maximum  temperature  excess  at  the  core  of  the  plasma,  r  is  the  radial  distance 
from  the  center  and  Tq  is  the  cut-off  radius  for  the  Gaussian  profile,  chosen  to 
be  tt/20  in  a  domain  of  < 

Energy  is  deposited  at  the  center  of  the  flow  domain  by  increasing  the  temper¬ 
ature  and  pressure  locally  at  constant  density.  The  initial  energy  deposited  A  E 
can  be  related  to  ATq  by  integrating, 

E  =  pG0cvAT4nr2dr,  (34) 

Jo 

which  yields 

E  —  ATon^poa^rl.  (35) 
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Figure  8:  Radial  profiles  of  density  at  t  —  0.014 r  0.034,  0.079,  0.141  and  0.218 
for  ATq  =  90  showing  shock  formation  and  propagation. 


Simulation  results  are  obtained  for  AT0  ranging  from  30  to  300.  Figure  7a  shows 
contours  of  pressure  at  t  =  1.2  for  a  maximum  temperature  excess  of  300  in  the 
core.  A  spherically  symmetric  shock  front  is  seen  to  propagate,  compressing 
the  flow  ahead  of  it,  while  the  flow  behind  it  expands.  Since  the  problem  is 
spherically  symmetric,  mean  values  of  the  flow  variables  are  computed  in  a 
radially  symmetric  manner,  and  plotted  as  a  function  of  radial  distance  from 
the  center  of  the  plasma. 

The  jump  in  flow  variables  across  the  shock  wave  will  be  incorrectly  predicted 
if  the  mesh  used  in  the  simulations  is  not  fine  enough  to  resolve  the  sharp 
gradients  at  the  shock  front.  A  grid  convergence  study  was  performed.  Figure 
7b  shows  the  radial  profile  of  density  at  the  instant  of  time  when  it  acquires  its 
maximum  value.  The  initial  temperature  ratio  in  the  core,  ATq  =  300.  The  first 
simulation  used  a  coarse  mesh  (Ax  —  0.0491).  The  mesh  size  was  then  refined 
successively  until  the  solution  did  not  change  significantly  for  two  successive 
mesh  refinements.  The  density  profiles  obtained  for  different  values  of  the  mesh 
spacing  are  shown  in  figure  7b.  The  solution  is  observed  to  be  grid  converged  for 
Ax  =  0.0196.  This  mesh  was  used  for  all  simulations  reported  in  this  section. 

3.1.1  Shock  formation  and  shock  propagation 

The  post-energy  deposition  process  can  be  divided  into  two  parts  -  shock  forma¬ 
tion  and  shock  propagation.  Figure  8  shows  radial  profiles  of  density  at  different 
instants  of  time.  Time  is  non-dimensionalized  with  respect  to  a  reference  time 
scale  fr  which  is  obtained  from  the  ratio  of  a  reference  length  scale  lr  and  refer¬ 
ence  velocity  scale  er>  For  the  spherical  idealization  tr  is  20/tt  times  the  initial 
radius  of  the  plasma  region,  while  cr  is  the  background  speed  of  sound.  There  is 
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Figure  9:  (a)  Radial  profiles  of  pressure  gradients  during  the  shock  formation 

process. - :  t  —  0.001, - :  t  —  0,005,  — —  :  f  =  0.013,  (b)  Radial 

profiles  of  pressure  gradients  during  the  shock  propagation  process. - :  t  — 

0.017, - :  t  =  0.023, - :  t  =  0.029. 


an  initial  period  of  time  during  which  the  intensity  of  the  shock  front  increases 
with  time,  indicating  formation  of  the  shock.  Once  shock  formation  is  com¬ 
plete,  the  shock  wave  propagates  radially  into  the  background  gas,  its  strength 
decreasing  as  a  result. 

Shock  formation  and  propagation  can  be  explained  by  the  large  pressure  gradi¬ 
ents  generated  by  energy  deposition.  Energy  is  deposited  in  the  form  of  internal 
energy  of  the  fluid  elements  by  increasing  the  temperature  at  constant  density. 
Since  the  initial  temperature  profile  is  a  Gaussian,  the  resulting  pressure  profile 
is  also  a  Gaussian.  The  gradients  in  pressure  thus  formed,  accelerate  the  flow 
radially,  through  the  conservation  of  momentum.  Thus  the  internal  energy  of 
the  fluid  elements  is  converted  into  kinetic  energy  through  the  pressure  gradients 
in  the  flow.  This  results  in  formation  of  a  shock  wave  whose  strength  initially 
increases.  The  pressure  gradients  then  drive  the  shock  wave  radially  into  the 
background  gas,  and  the  strength  of  the  shock  decreases. 

Figure  9a  shows  radial  profiles  of  pressure  gradients  at  different  instants  dur¬ 
ing  shock  formation.  Note  that  the  gradients  in  pressure  are  initially  negative, 
indicating  that  the  pressure  is  highest  at  the  center  of  the  plasma.  At  later 
instants,  the  pressure  gradients  acquire  positive  and  negative  values,  indicating 
that  a  shock  front  has  begun  to  form,  and  the  peak  value  of  pressure  is  no  longer 
at  the  center  of  the  domain.  Figure  9b  shows  radial  profiles  of  the  pressure  gra¬ 
dients  as  the  shock  propagates.  Pressure  gradients  at  three  subsequent  instants 
of  time  are  shown.  Note  that  both  the  positive  and  negative  values  of  pressure 
gradients  decrease  in  time  due  to  the  radial  spreading  of  the  shock  wave. 

Figure  10a  shows  the  variation  of  shock  Mach  number  in  time  for  ATq  =  300, 
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Figure  10:  (a)  Variation  of  shock  Mach  number  with  time  - :  ATq  —  300, 

—  —  :  A7o  =  210  , - :  A  Jo  =  90,  (b)  Variation  of  maximum  shock 

Mach  number  with  Alb,  a  :  simulation  results,  ■  :  Experimental  data  (Van 
et  al  2003). 


210  and  90.  The  initial  steep  increase  in  the  shock  Mach  number  depicts  shock 
formation.  The  later  gradual  decrease  in  shock  Mach  number  depicts  the  radial 
propagation  of  the  shock.  The  shock  Mach  number  attains  a  maximum  value 
at  the  end  of  shock  formation.  The  values  of  shock  Mach  number  increase  with 
increase  in  energy  deposited. 

Figure  10b  shows  the  maximum  shock  Mach  number  for  different  values  of  ATq. 
The  shock  Mach  number  can  be  computed  using  the  relation 

U,  =  +  ({lilt/)3 +C*3)1'3,  (36) 

for  the  shock  velocity  U9.  Here,  U  is  the  velocity  behind  the  shock  front  and 
Coo  is  the  speed  of  sound  ahead  of  it  The  maximum  shock  Mach  number 
is  plotted  for  ATo  ranging  from  30  to  300.  As  expected,  M9  increases  with 
energy  deposited  in  the  How.  A  power  law  fit  is  obtained  and  extrapolated 
to  higher  temperatures.  The  prediction  of  maximum  shock  Mach  number  for 
higher  temperatures  matches  well  with  data  from  experiments  at  AT0  —  570 
(Yan  et  al  2003). 

Figure  11a  shows  shock  formation  time  as  a  function  of  AT0.  Figure  lib  shows 
shock  propagation  time  for  different  values  of  A 7b.  Here,  shock  propagation 
time  is  the  time  required  for  the  shock  front  to  reach  a  point  where  r/r0  — 
40/tt,  r0  being  the  initial  size  of  the  plasma.  Both  shock  formation  and  shock 
propagation  time  decrease  with  increase  in  ATq,  Comparison  of  figures  11a  and 
lib  shows  that  for  a  given  AT0  the  shock  propagation  time  is  2  to  3  orders  of 
magnitude  larger  than  the  corresponding  shock  formation  time.  This  can  also 
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Figure  11:  (a)  Shock  formation  and  (b)  Shock  propagation  time  for  different 

AT0. 


be  observed  in  figure  10.  In  fact,  for  A  To  =  300T  the  shock  formation  time 
becomes  comparable  to  the  time  scale  in  which  the  plasma  generation  occurs 
and  so  the  Navier- Stokes  model  may  not  accurately  represent  some  of  the  flow 
features. 

3*1*2  Comparison  to  strong  shock  solution 

Based  on  dimensional  arguments,  Taylor  (1950)  and  Sedov  (1959)  independently 
proposed  a  similarity  solution  for  the  point  blast  problem.  The  process  was 
idealized  as  a  sudden  release  of  energy  E  at  a  point.  Further,  it  was  assumed 
that  the  resulting  shock  wave  is  so  strong  that  the  pressure  and  sound  speed  of 
the  background  medium  are  negligible  compared  to  the  pressure  and  velocity  of 
the  shock  wave.  Based  on  these  assumptions,  the  only  dimensional  parameters 
involved  in  this  problem  are  the  amount  of  energy  deposited  Er  and  the  density 
of  the  background  medium  The  strong  shock  jump  conditions  (Taylor  1950, 
Sedov  1959)  apply;  he. 

2  -j-  1  2  <i 

u  ~  '  p  ~  poa  1  p  =  ThTT Pao Us  (■ 37) 

where  u,  p  and  p  are  the  velocity,  density  and  pressure  immediately  behind  the 
shock  wave,  and  Ua  is  the  shock  velocity.  Then  based  on  dimensional  arguments, 
the  shock  radius  Rs  at  any  time  t  is  given  by 

/  E*  \ 

R,(t)  =  fe(— J  *2/S  (38) 
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Figure  12:  (a)  Variation  of  shock  radius  with  time,  (b)  Variation  of  /^/(ATq)1/5 

with  time,  - :  A7b  —  300,  •  :  ATo  =  240,  « :  A2o  -  180, - :  ATq  ~ 

30. 


where  fc  is  a  dimensionless  number.  The  pressure  and  velocity  immediately 
behind  the  shock  front  can  be  obtained  from  the  jump  conditions 


P 


u 


8  (JL\2'Sre/ s 

257  +  1  \PooJ 

57  +  1  \PcoJ 


(39) 

(40) 


Simulation  results  for  variation  of  shock  radii  in  time  for  AT0  —  300,  240, 
180  and  30  are  shown  in  figure  1 2a.  Note  that  from  equation  35, 


E  ATq,  (41) 

Equations  41  and  38  suggest  that  if  the  strong  shock  solution  were  to  apply,  then 
all  the  curves  in  figure  12a  should  collapse  into  one  single  curve  when  R9/ATq1^ 
is  plotted  against  time.  Figure  12b  shows  the  variation  of  R9/ATq1^  with  time 
for  A  Jo  =  300,  240,  180  and  30.  It  is  clear  that  for  an  initial  time  all  these  curves 
collapse  over  one  another,  but  with  increasing  time  they  start  to  separate.  The 
lowest  temperature  ratio  is  the  first  to  deviate.  The  curves  would  be  expected 
to  collapse  only  as  long  as  the  strong  shock  conditions  are  applicable.  Smaller 
the  energy  deposited,  sooner  the  shock  strength  weakens,  and  sooner  the  curve 
deviates  from  the  strong  shock  limit. 

Note  that  deviation  from  the  strong  shock  limit  leads  to  larger  values  of  Rs/ATq1/^ 
than  would  be  expected  from  the  strong  shock  solution.  This  can  be  explained 
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Figure  13:  (a)  Variation  of  velocity  at  the  shock  front  against  time,  (b)  Variation 
of  w/tATo)1^5  against  time,  - :  A T0  =  300,  •  :  AT0  =  210,  ■ :  A7b  —  90. 


as  follows.  As  the  shock  propagates  radially,  its  strength  decreases,  and  even¬ 
tually  it  becomes  an  acoustic  wave  traveling  at  the  speed  of  sound.  Hence  the 
other  relevant  limit  for  the  time  variation  of  shock  radius  is  the  acoustic  limit 
which  the  shock  wave  eventually  attains,  regardless  of  the  energy  deposited  in 
the  flow.  The  time  required  to  reach  this  limit,  would  however  increase  as  the 
energy  deposited  increases.  For  the  acoustic  limit, 


(42) 


where  Coq  is  the  background  speed  of  sound.  Thus  for  the  curves  in  figure  12b, 
the  slope  in  the  acoustic  limit  would  be  Coc/ATb1^5.  So,  for  smaller  energy 
deposited  in  the  flow  the  slope  of  the  curve  is  larger  in  the  acoustic  limit,  and 
hence  the  deviation  from  the  strong  shock  limit  is  towards  larger  values  of 

RJAT0i/s. 

Again,  combining  equations  41  and  40,  it  may  be  reasoned  that  in  the  strong 
shock  limit,  the  curves  of  u/ATq 3j/5  against  time  must  collapse  over  one  another 
for  different  values  of  ATq,  Figure  13a  shows  velocity  immediately  behind  the 
shock  front  against  time  for  ATb  =  300,  210  and  90,  Figure  13b  shows  uj  AT0I/5 
against  time  for  A7b  —  300,  210  and  90.  Collapse  with  the  strong  shock  scal¬ 
ing  is  reasonable  at  initial  times  when  the  strong  shock  conditions  are  valid. 
However,  at  later  instants  of  time,  these  curves  deviate  from  the  strong  shock 
solution.  This  behavior  can  again  be  explained  in  terms  of  deviation  from  the 
strong  shock  limit  towards  the  acoustic  limit. 
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Figure  14:  (a)  Schematic  to  illustrate  the  conservation  of  mass  within  a  con¬ 
trol  volume  surrounding  the  shock  wave,  (b)  Minimum  density  in  the  core  for 
different  amounts  of  energy  deposited  in  the  flow. 


3.1.3  Behavior  of  the  core 

The  strong  shock  solution  solution  predicts  the  density  at  the  core  to  be  zero. 
This  would  mean  that  the  core  has  expanded  to  the  extent  that  there  is  no 
mass  left  in  the  core.  However  for  weaker  conditions,  as  the  density  at  the 
shock  front  keeps  increasing  during  shock  formation,  the  density  in  the  core 
keeps  decreasing.  Once  the  shock  formation  is  complete,  the  shock  propagates 
down  the  domain  with  decreasing  intensity.  Hence  density  at  the  shock  front 
starts  decreasing  and  the  density  at  the  core  simultaneously  starts  increasing 
(Figure  8). 

This  behavior  can  be  explained  by  conservation  of  mass  behind  the  shock  front. 
Figure  14a  shows  a  schematic  for  density  variation  behind  the  shock  front.  Con¬ 
sider  a  spherical  control  volume  of  radius  R0  around  the  shock- wave  whose 
radius  R$  <  Rq.  The  integral  form  of  the  continuity  equation  is  given  by 


f  \dp  dpuj 

Jv  [eft  dxj  m 


dV  =  0 


(43) 


where  V  is  the  control  volume  under  consideration.  Since  there  is  no  mass  flux 
in  or  out  of  the  control  volume,  the  second  term  in  equation  43  becomes  zero. 
Also  assuming  that  the  control  volume  is  targe,  the  shock  front  is  confined  within 
this  volume  for  all  time  and  so  the  control  volume  does  not  change  with  time. 
Then  the  time  derivative  in  the  first  term  of  equation  43  can  be  taken  outside 
the  volume  integral  to  obtain 


=  m0 


(44) 
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Figure  15:  (a)  Radial  profiles  of  pur  at  different  instants  of  time  showing  radial 

collapse  of  the  flow  field,  - :  t  =  0.051,  - - :  t  -  0.186, - :  t  « 

0,422,  (b)  Radial  location  of  the  spherical  stagnation  front  in  time. 
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Figure  16:  (a)  Radial  profiles  of  internal  energy  per  unit  volume  at  different 
instants  of  time,  (b)  Radial  profiles  of  kinetic  energy  per  unit  volume  at  dif¬ 
ferent  instants  of  time,  -  :  t  —  0.005,  - - :  t  —  0.011, - -  :  t  — 

0.051 1 .  :  t  =  0.186. 
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Figure  17:  (a)  Axial  variation  of  the  initial  temperature  profile,  (b)  Contour 
plot  of  initial  temperature  in  a  plane  through  the  plasma  axis. 


where  m0  is  a  constant  and  is  equal  to  the  initial  mass  within  the  control  volume. 
Thus  total  mass  within  the  control  volume  is  conserved  at  all  times. 

At  the  shock  front,  mass  accumulates  due  to  compression.  This  mass  comes 
from  the  core,  and  hence  the  core  experiences  an  expansion  leading  to  decrease 
in  density.  During  shock  formation,  mass  accumulation  at  the  shock  front  in¬ 
creases  and  so  density  in  the  core  decreases.  The  reverse  occurs  during  shock 
propagation.  At  any  instant  of  time  the  accumulated  mass  at  the  shock  front 
indicated  by  the  shaded  area  A  has  to  balance  the  deficit  of  mass  in  the  core 
indicated  by  the  shaded  area  B. 

Figure  14b  shows  the  minimum  value  of  density  at  the  core  as  a  function  of 
energy  deposited.  The  minimum  value  of  density  is  observed  to  decrease  with 
increase  in  energy  deposited  in  the  Row.  The  trend  suggests  that  in  the  strong 
shock  limit,  the  density  in  the  core  will  tend  to  be  zero.  However  with  large 
decrease  in  density,  note  that  the  continuum  model  used  for  the  Navier-Stokes 
equations  becomes  less  valid. 

Note  that  as  the  shock  Front  propagates  out  radially,  the  flow  behind  the  shock 
front  experiences  a  radial  collapse.  Figure  1 5a  shows  radial  profiles  of  the  con¬ 
vection  flux,  puri  at  different  instants.  The  initial  temperature  ratio  in  the  core, 
ATb  ”  90.  During  formation  of  the  shock  wave,  the  flow  expands  outwards  and 
no  reverse  flow  is  observed.  However  as  the  shock  wave  propagates  into  the  back¬ 
ground  medium,  the  outward  radial  motion  near  the  shock  front  is  followed  by  a 
reverse  radial  flow  in  the  core.  Hence  a  spherical  stagnation  front  is  formed  that 
separates  the  flow  moving  radially  inwards  from  that  moving  radially  outwards. 
The  location  of  the  stagnation  front  depends  on  the  location  of  the  shock  front 
and  moves  outward  in  time.  Figure  15b  shows  radial  location  of  the  spherical 
stagnation  front,  Rspi  as  a  function  of  time.  The  initial  flat  portion  of  the  curve 
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Figure  18:  {a)  Comparison  of  shock  radius  normal  to  the  axis  of  the  plasma,  (b) 

Density  profiles  normal  to  the  axis  of  the  plasma  at  10  ps  - :  simulation 

results, - :  experiment  {Adelgren  et  al  2001),  - - :  simulations  of 

Kandala  and  Candler  (2003). 


represents  formation  of  the  shock  wave  when  the  flow  expands  radially  outward 
only.  As  the  reverse  flow  sets  in  the  radial  location  of  the  stagnation  front  is 
observed  to  vary  almost  linearly  with  time. 

Figure  16a  shows  radial  profiles  of  internal  energy  per  unit  volume,  pe}  at  dif¬ 
ferent  instants  of  time.  Radial  profiles  of  kinetic  energy  per  unit  volume,  pk , 
are  shown  in  figure  1 6b.  The  initial  temperature  ratio  in  the  core,  AT0  =  90. 
During  formation  of  the  shock  wave,  denoted  by  the  first  two  profiles  in  figure 
16t  internal  energy  of  the  fluid  elements  is  observed  to  get  converted  into  ki¬ 
netic  energy  leading  to  high  velocities  at  the  shock  front.  However  as  the  shock 
propagates  radially,  the  velocities  at  the  shock  front  decrease  and  the  kinetic 
energy  is  converted  back  to  internal  energy  as  observed  in  the  last  two  profiles 
in  figure  16,  Note  that  within  the  control  volume  shown  in  figure  14a,  the  total 
energy  is  always  conserved. 

3.2  Tear-drop  shaped  idealization 

This  section  contains  results  from  simulations  for  the  tear-drop  shaped  ideal¬ 
ization  of  the  plasma.  Experiments  in  laser  induced  breakdown  (Adelgren  et  al 
2001,  Yan  et  al  2003)  show  that  the  plasma  is  initially  tear-drop  shaped.  Figure 
17a  shows  the  axial  temperature  distribution  that  is  used  to  model  the  initial 
temperature  profile  for  the  tear-drop  idealization  of  the  plasma.  Laser  energy 
deposition  is  assumed  to  be  symmetric  about  the  laser  axis.  The  temperature 
profile  normal  to  the  plasma  axis  is  assumed  to  be  a  Gaussian.  Figure  17b  shows 
contours  of  temperature  in  a  plane  passing  through  the  axis  of  the  plasma.  The 
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Figure  19;  (a)  Pressure  profiles  normal  to  the  axis  of  the  plasma  at  10  /is, 

(b)  velocity  profiles  normal  to  the  axis  of  the  plasma  at  10  fis  - ;  simulation 

results, - :  experiment  (Adelgren  et  al  2001),  - - :  simulations  of 

Kamlala  and  Candler  (2003). 


initial  temperature  profile  is  similar  to  that  observed  in  simulations  of  Kandala 
and  Candler  (2003). 

3.2.1  Comparison  to  existing  results 

Figure  1 8a  compares  the  shock  radius  normal  to  the  plasma  axis  with  experiment 
(Adelgren  et  al  2001)  and  simulation  (Kandala  and  Candler  2003).  Shock  radius 
is  defined  by  the  outer  edge  of  the  pressure  blast  wave.  Here  the  reference  length 
scale  lr  is  obtained  from  the  initial  length  of  the  plasma  region  to  be  4.4075  mm 
and  the  reference  velocity  scale,  cr,  is  given  by  the  background  speed  of  sound. 
The  reference  time  scale,  tTy  is  thus  obtained  to  be  12.695  /is.  The  shock  radius 
is  seen  to  be  slightly  under-predicted.  Figures  18b,  19a  and  19b  show  pressure, 
density  and  velocity  profiles  normal  to  the  axis  of  the  plasma.  The  profiles  are 
observed  to  peak  slightly  behind  the  peaks  obtained  in  experiments.  However, 
the  peak  values  of  density,  pressure  and  velocity  are  well  predicted.  Also  the 
behavior  of  density  in  the  core  is  captured  well.  The  reason  for  the  differences 
may  be  attributed  to  the  lack  of  accurate  information  about  the  variation  of 
the  initial  temperature  profile  normal  to  the  axis  of  the  plasma  due  to  which 
the  temperature  has  been  assumed  to  vary  as  a  Gaussian  in  this  direction. 
The  reasonable  agreement  with  existing  data  suggests  that  for  the  purpose  of 
studying  the  properties  of  the  blast  wave  and  its  effect  on  the  surrounding  fluid, 
the  Navier-Stokes  idealization  might  be  adequate. 
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Figure  20:  Pressure  contours  at  t  =  0,02  and  0,03  show  formation  of  a  distinct 
shock  wave  normal  to  the  axis  of  the  plasma. 
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Figure  21:  Pressure  contours  at  t  =  0.36  and  1.5  show  that  the  blast  wave  is 
initially  tear-drop  shaped  but  become  spherical  in  time. 
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Figure  22:  (a)  Schematic  showing  propagation  of  a  tear-drop  shaped  blast 
wave,  (b)  Generation  of  vorticity  near  leading  and  trailing  edges  of  the  plasma, 

-  :  vorticity  jump  across  the  shock  wave,  — —  :  shock  strength, 

- :  curvature  of  the  shock  front,  »  :  tangential  velocity  measured  ahead 

of  the  shock. 


3,2,2  Shock  formation  and  propagation 

As  done  for  the  spherical  idealization,  laser  energy  is  added  into  the  flow  by 
increasing  the  temperature  at  constant  density.  Thus  sharp  gradients  in  pressure 
are  generated  in  the  flow.  Again  the  entire  process  can  be  divided  into  the 
formation  and  propagation  of  a  shock  wave.  Formation  of  the  shock  wave  is 
shown  in  figure  20.  Contours  of  pressure  are  shown  at  two  successive  instants 
of  time.  Due  to  the  tear-drop  shape  of  the  plasma,  pressure  gradients  normal 
to  the  plasma  axis  are  larger  compared  to  the  gradients  along  the  plasma  axis. 
Thus  fluid  elements  in  a  region  normal  to  the  axis  begin  to  gain  kinetic  energy 
and  a  shock  wave  begins  to  form,  and  gradually  steepen  in  intensity.  The 
shock  wave  then  propagates  into  the  background  gas  with  decreasing  intensity. 
Figure  21  shows  contours  of  pressure  at  two  different  instants  of  time  during 
propagation  of  the  shock  wave.  The  blast  wave  is  initially  tear-drop  shaped  as 
observed  in  figure  21a.  The  stronger  pressure  gradients  normal  to  the  plasma 
axis  leads  to  larger  acceleration  of  fluid  elements  in  this  direction.  The  effect  is 
to  push  the  blast  wave  faster  normal  to  the  plasma  axis  until  pressure  gradients 
in  all  directions  become  uniform  and  a  nearly  spherical  blast  wave  propagates 
into  the  domain.  Figure  21b  shows  that  as  the  blast  wave  propagates  to  a  radius 
2  to  3  times  the  initial  size  of  the  plasma  it  becomes  more  and  more  spherical 
in  shape. 
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Figure  23:  Contours  of  $  component  of  vorticity  at  t  =  0.02  and  0.08, 


:  positive  vorticity, - :  negative  vorticity. 


3,2,3  Short  time  behavior 


An  important  feature  of  the  tear-drop  shaped  idealization  is  generation  of  baro- 
dinic  vorticity  near  the  leading  and  trailing  edges  of  the  plasma  for  times  less 
than  5  ps.  Shock  waves  with  varying  strength  and  curvature  generate  vorticity 
in  the  flow  (TYuesdell  1952,  Lighthill  1957,  Hayes  1957,  Berndt  1966,  Kevlahan 
1997),  A  general  expression  for  vorticity  jump  across  an  unsteady  curved  shock 
propagating  into  a  two-dimensional  non-uniform  flow  was  derived  by  Kevlahan 
(1997).  However  a  simpler  expression  for  vorticity  jump  suggested  by  TYuesdell 
(1952)  has  been  used  to  explain  generation  of  vorticity  near  the  leading  and 
trailing  edges  of  the  plasma.  TYuesdell  considered  vorticity  jump  across  a  two 
dimensional  steady  shock  in  uniform  flow.  The  jump  in  vorticity  was  given  by 


Suj  = 


1  +/* 


nUu 


(45) 


where  k  is  the  curvature  of  the  shock  wave,  Ut  is  the  velocity  tangential  to  the 
shock  in  the  frame  of  reference  of  the  shock  wave,  and  p  is  the  shock  strength 
given  by 


n  —  — —  ii  {46} 

poo 

where  ps  and  poo  are  the  densities  at  the  shock  front,  and  ahead  of  it  respectively. 
Note  that  vorticity  has  been  considered  to  be  positive,  out  of  the  plane  of  the 
paper. 

Figure  22a  shows  a  schematic  for  propagation  of  a  tear-drop  shaped  shock  wave. 
Note  that  both  curvature  and  strength  of  the  shock  wave  varies  along  its  surface. 
Figure  22b  shows  the  variation  of  shock  strength,  shock  curvature,  tangential 
velocity  and  jump  in  vorticity  as  a  Function  of  z  at  t  =  0.08.  The  curvature  of 
the  shock  wave  is  zero  near  the  center  and  obtains  large  negative  values  near  the 
edges  of  the  plasma.  The  shock  strength  on  the  other  hand  is  maximum  near  the 
center  and  decreases  towards  the  edges.  The  tangential  velocity  is  negative  and 
positive  near  the  leading  and  trailing  edges.  The  jump  in  vorticity  as  obtained 
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Figure  24:  (a)  Baroclinic  term  as  a  function  of  z  at  r  =  0.172  and  at  t  ~  0.02, 
(b)  Bulk  dilatation  term  as  a  Function  of  z  at  r  =  0.172  and  at  t  =  0.02,  LE  and 
TE  denote  leading  and  trailing  edges  of  the  plasma. 


From  equation  45,  is  hence  found  to  be  positive  and  negative  near  the  leading 
and  trailing  edges  of  the  plasma,  respectively. 

The  vorticity  equation  for  a  compressible,  viscous  flow  with  variable  fluid  prop¬ 
erties  can  be  written  as 


dd) 

dt 


(u  ■  V)l3  +  (w  ■  V)u 

convection  vortex  stretching/tilting 

Vp  x  Vp 


-  d?(V  -  u) 

bulk  dilatation 

2-(Vx(iv.T))(47) 
/Ce  p 


baroclinic 


Since  baroclinic  vorticity  is  the  only  source  term  in  equation  47,  it  is  the  only 
term  responsible  for  initial  generation  of  vorticity  in  the  flow.  Note  that  the 
flow  is  axisymmetric.  Defining  a  cylindrical  coordinate  system  along  rt  $  and  z 
(figure  22a),  any  gradients  in  the  flow  will  occur  only  along  f  and  z  directions. 
Thus  the  only  possible  component  of  baroclinic  vorticity  in  the  flow  will  be  along 

l 

Figure  23  shows  contours  of  the  0  component  of  vorticity  at  subsequent  instants 
of  time.  The  contours  are  shown  in  a  plane  passing  through  the  axis  of  the 
plasma.  Positive  and  negative  vorticity  is  observed  near  the  leading  and  trailing 
edges  of  the  plasma.  Vorticity  generation  is  stronger  near  the  leading  edge.  The 
magnitude  of  vorticity  is  observed  to  decrease  in  time,  as  it  spreads  out  further 
from  the  core  of  the  plasma. 

Budgets  were  computed  for  the  different  terms  on  the  right  hand  side  of  the 
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(a) 


(b) 


Figure  25:  (a)  Vorticity  budget  at  a  point  near  the  leading  edge  of  the  plasma 
(z  =  3.5,  r  =  0.172),  (b)  Vorticity  budget  at  point  near  the  trailing  edge  (z  — 

2.8,  r  =  0*172).  - :  convection,  — —  :  vortex  stretching/ til  ting  , 

- -  :  bulk  dilatation,  - :  baroclinie,  . .  :  viscous. 


vorticity  equation  and  their  behavior  in  time  was  observed.  Figure  24a  shows  the 
baroclinie  source  term  as  a  function  of  z  at  r  =  0.172  and  at  t  —  0.02.  Positive 
and  negative  regions  of  baroclinie  vorticity  are  observed  near  the  leading  and 
trailing  edges  of  the  plasma  respectively*  The  region  with  positive  vorticity  is 
stronger  compared  to  that  with  negative  vorticity.  This  is  because  density  is 
lower  near  the  leading  edge  of  the  plasma.  Also  the  curvature  of  the  blast  wave 
changes  rapidly  near  the  leading  edge  causing  maximum  misalignment  between 
the  pressure  and  density  gradients.  Note  that  from  symmetry  considerations, 
the  alignment  of  Vp  and  Vp  near  the  trailing  edge  is  just  the  opposite  to  that 
near  the  leading  edge.  So  the  baroclinie  source  term  is  opposite  in  sign  near  the 
leading  and  trailing  edges  of  the  plasma.  Figure  24b  shows  the  bulk  dilatation 
term  as  a  function  of  z  at  r  —  0.172  at  the  same  instant  of  time.  Since  there  is 
expansion  behind  the  blast  wave  (V  u)  is  positive.  Since  vorticity  is  positive 
and  negative  near  the  leading  and  trailing  edges  of  the  plasma  (figure  23),  the 
bulk  dilatation  term  is  negative  and  positive  near  the  leading  and  trailing  edges 
respectively. 

Figure  25  shows  time  variation  of  the  different  terms  on  the  right  hand  side 
of  the  vorticity  equation.  Near  the  leading  edge  of  the  plasma  {figure  25a), 
the  generation  of  baroclinie  vorticity  is  positive.  This  generation  of  vorticity  is 
suppressed  in  time  through  the  bulk  dilatation  term  which  is  negative  near  the 
leading  edge  of  the  plasma.  The  reverse  occurs  near  the  trailing  edge  (figure 
25b)  where  the  alignment  of  net  pressure  and  density  gradients  is  such  that  the 
baroclinie  vorticity  generated  is  negative.  The  bulk  dilatation  term  is  positive 
near  the  trailing  edge  and  in  time  counter  balances  the  baroclinie  generation  of 
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Figure  26:  Temperature  contours  at  t  =  18  ps,  30  /is  and  60  fis  show  breaking 
and  roll  up  of  the  temperature  field. 


vorticity. 

Vorticity  is  observed  to  be  convected  in  the  direction  of  propagation  of  the  tear¬ 
drop  shaped  blast  wave.  Thus  the  spread  of  vorticity  away  from  the  plasma 
axis  (figure  23)  can  be  attributed  to  the  combined  effects  of  vorticity  convection 
and  propagation  of  the  blast  wave  resulting  in  vorticity  generation  at  locations 
further  from  the  axis  of  the  plasma  The  vortex  stretching/ tilting  term  is  found 
to  be  positive  and  negative  near  the  leading  and  trailing  edges  respectively.  The 
viscous  term  is  2  to  3  orders  of  magnitude  smaller  than  any  of  the  other  terms 
and  so  does  not  have  a  significant  contribution  in  the  evolution  of  vorticity  in 
the  flow.  So,  vorticity  generated  by  baroclinic  means  due  to  interaction  of  the 
curved  blast  wave  with  the  flow  upstream,  is  convected  radially  outward,  and 
suppressed  by  bulk  expansion  in  the  core  of  the  plasma. 

3.2.4  Long  time  behavior 

This  section  shows  evolution  of  the  core  of  the  plasma  at  times  longer  than  5  fts. 
Experiments  in  laser-induced  break  down  (  Adelgren  2001,  Yan  2003)  show  that 
the  core  of  the  plasma  eventually  roll  up  generating  a  toroidal  vortex  ring.  The 
time  scale  for  this  roll  up  depends  on  the  amount  of  laser  energy  deposited  in  the 
flow.  Experimental  temperature  measurements  indicate  that  as  a  consequence 
of  the  roll  up,  the  temperature  field  is  advected  such  that  the  peak  temperature 
shifts  from  the  plasma  axis  to  the  center  of  the  generated  vortex  ring. 

Figure  26  shows  long  time  evolution  of  the  temperature  field  obtained  in  the 
simulations.  Note  that  the  temperature  field  starts  breaking  around  t  —  18  fis 
(figure  26a),  The  break  up  starts  from  the  region  of  maximum  temperature 
along  the  axis  of  the  plasma.  This  propagates  along  the  axis  of  the  plasma  from 
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Figure  27:  Velocity  streamlines  at  t  —  4  jjs,  12  /is  and  30  [is. 


Figure  28:  Velocity  streamlines  plotted  over  contours  of  vortidty  magnitude  at 
30 
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Figure  29:  Contours  of  density,  pressure  and  temperature  at  a  given  instant  of 
time  shows  interaction  of  the  spherical  plasma  with  a  turbulent  background. 


right  to  left  resulting  in  formation  of  an  ax isym metric  temperature  lobe  around 
30  {is  (figure  26b).  In  time  this  temperature  lobe  moves  further  away  from  the 
plasma  axis  and  finally  rolls  up  to  form  a  toroidal  vortex  ring  as  shown  in  the 
contours  of  temperature  at  60  {is  (figure  26c),  Note  that  during  this  process, 
the  maximum  temperature  is  advecteci  from  the  plasma  axis  to  the  center  of  the 
vortex  ring. 

The  breaking  and  roll  up  of  the  temperature  field  can  be  explained  by  evolution 
of  the  velocity  field  in  the  core  of  the  plasma.  Figure  27  shows  time  evolution 
of  the  velocity  field.  Initially  the  core  of  the  plasma  expands.  This  is  observed 
from  velocity  streamlines  moving  away  from  the  plasma  core  (figure  27a),  As  the 
blast  wave  propagates  away  from  the  core,  it  can  be  shown,  from  conservation  of 
mass  that  part  of  the  fluid  from  the  shock  front  moves  back  towards  the  plasma 
core.  Thus  a  separating  streamline  is  formed,  ahead  of  which  the  fluid  elements 
move  away  from  the  plasma  core,  but  behind  which  a  reverse  flow  is  observed 
(figure  27b).  The  reverse  flow  is  strongest  along  the  plasma  axis  from  right  to 
left.  Thus  two  jets  of  different  velocity  impinge  head  on  along  the  plasma  axis, 
yielding  a  complicated  vorticity  field  (figure  27c).  The  fact  that  fluid  rushes  in 
from  right  to  left  along  the  plasma  axis  and  then  rolls  over  in  a  direction  normal 
to  the  plasma  axis  leads  to  breaking  and  roll  up  of  the  temperature  field. 

Note  that  the  actual  vorticity  field  is  much  more  complicated  than  just  a  single 
toroidal  vortex  ring.  However  the  magnitude  of  vorticity  is  much  stronger  in  the 
vortex  ring  closest  to  the  plasma  axis.  This  can  be  shown  by  plotting  velocity 
streamlines  over  contours  of  vorticity  magnitude  (figure  28),  Thus  only  a  single 
toroidal  vortex  ring  is  observed  in  the  experimental  flow  visualization. 
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Figure  30;  Statistics  for  the  radial  component  of  velocity  and  transverse  com¬ 
ponent  of  vorticity  at  three  different  instants  of  time. 

3.3  Spherical  plasma  interacting  with  isotropic  turbu¬ 
lence 

Direct  numerical  simulation  of  isotropic  turbulence  was  performed  under  condi¬ 
tions  corresponding  to  past  DNS  by  Blaisdell  et  ah  (1991).  The  initial  velocity 
fluctuations  are  isotropic  and  divergence-free,  while  initial  fluctuations  in  pres¬ 
sure,  temperature  and  density  are  assumed  to  be  zero.  The  initial  velocity  fluc¬ 
tuations  are  generated  using  Rogallo's  (1981)  method.  The  Fourier  coefficients 
of  the  initial  velocity  fields  are  given  by 


akk2  +  0k\ k$  0^2^  -  akki  0(k\2  -F  fca2)1^3 

u  =  k(k  TTkTyji6' +  fc<fc>a  +  *2!)^ej  k  C3 


(48) 


where  ^  ^ 

=  ei9lcosW-  e'e,cos(<£),  (49) 

Here  $u  O2  and  <t>  are  random  numbers  from  0  to  2it,  and  e*  denote  the  unit 
vectors  along  the  three  coordinate  directions. 

The  initial  energy  spectrum  used  in  the  simulations  is  a  power  four  spectrum 
peaking  at  kQ  =  5,  The  initial  flow  field  thus  generated  has  an  fluctuation  Mach 
number  Mt  -  0,3  and  turbulent  Reynolds  number  Rex  =  50.  The  turbulent 
field  is  allowed  to  decay  for  some  time  after  which  the  velocity  derivative  skew¬ 
ness  attains  a  steady  value  in  the  range  of  “0.3  to  -0,4.  The  plasma  is  then 
introduced  by  increasing  the  temperature  and  pressure  at  constant  density. 
Figure  29  shows  contours  of  density,  pressure  and  temperature  at  a  later  instant 
of  time.  A  spherically  symmetric  blast  wave  is  observed  to  propagate  through 
the  background  turbulence  compressing  the  flow  in  a  spherically  symmetric 
manner.  There  is  significant  expansion  in  the  plasma  core.  As  the  blast  wave 
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Figure  31:  Instantaneous  isosurfaces  of  velocity  magnitude  showing  propagation 
of  a  blast  wave  through  isotropic  turbulence. 


propagates  significantly  into  the  background,  it  is  observed  to  get  distorted 
rapidly  due  to  interaction  with  the  turbulence. 

Since  the  problem  has  radial  symmetry,  statistics  were  obtained  in  a  radially 
symmetric  manner.  Figure  30  shows  statistics  for  the  radial  component  of  ve¬ 
locity  and  the  transverse  component  of  vortidty  at  three  different  instants  of 
time.  Amplification  of  turbulence  levels  are  observed  in  the  vicinity  of  the  blast 
wave.  The  local  compression  across  the  blast  wave  causes  the  turbulence  levels 
to  amplify.  Also  since  the  blast  wave  is  radially  symmetric  the  amplification 
is  most  pronounced  in  the  radial  component  of  the  velocity  and  hence  in  the 
transverse  component  of  the  vorticity. 

3.4  Tear-drop  shaped  plasma  interacting  with  isotropic 
turbulence 

Laser  energy  is  added  to  this  flow  by  locally  increasing  the  temperature  and 
pressure  at  constant  density.  Figure  31  shows  a  schematic  for  the  problem. 
Isosurfaces  of  velocity  magnitude  has  been  shown  at  a  given  instant  of  time.  A 
tear-drop  shaped  blast  wave  is  observed  to  propagate  through  the  background 
turbulence.  Interaction  of  this  blast  wave  with  the  background  turbulence  and 
the  effect  of  expansion  in  the  plasma  core  has  been  described  in  this  section. 
Figure  32  shows  propagation  of  the  blast  wave  through  background  turbulence. 
Both  side  and  sectional  views  have  been  shown  for  contours  of  density  plotted 
at  a  giiTen  instant  of  time.  A  very  strong  blast  wave  is  observed  to  propagate 
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Figure  32:  Contours  of  density  at  t  =  275  ps  showing  propagation  of  a  strong 
blast  wave  through  the  background  turbulence. 


z  z  z 

Figure  33:  Statistics  for  turbulent  kinetic  energy,  vorticity  and  divergence  at 
t  =  L9  its. 


into  the  background  compressing  the  turbulence,  leaving  the  flow  in  the  core 
to  expand.  The  blast  wave  is  initially  strong  and  hence  retains  its  tear  drop 
shape.  But  the  effect  of  the  background  turbulence  is  to  dampen  the  blast 
wave  intensity  and  soon  it  becomes  weak  enough  to  interact  heavily  with  the 
turbulence  and  its  shape  gets  distorted. 

Compression  across  the  blast  wave  and  expansion  in  the  plasma  core  affect 
turbulence  levels  in  these  regions.  Note  that  the  turbulence  is  statistically  ho¬ 
mogeneous  in  axisym  metric  planes  with  respect  to  the  plasma  axis.  Statistics 
were  therefore  computed  in  z  —  r  planes.  Figure  33  shows  statistics  for  turbulent 
kinetic  energy,  vorticity  and  divergence  at  t  —  1.9  iis.  Note  that  in  the  vicinity 
of  the  blast  wave  turbulence  levels  get  amplified  due  to  presence  of  compression 
there.  This  trend  is  supported  by  linear  analysis  of  turbulence  interacting  with 
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Figure  34:  Statistics  for  turbulent  kinetic  energy  and  divergence  at  t  =  1.9  fis, 
2.25  ps  and  2.5  fis. 

a  shock  wave  {Mahesh  et  al.  1995,1997).  Also  the  expansion  in  the  plasma  core 
leads  to  suppression  of  turbulence  levels  there.  Increase  in  viscosity  due  to  high 
temperatures  in  the  core  also  add  to  suppression  of  turbulence  levels. 

The  extent  of  turbulence  amplification  across  a  shock  wave  would  depend  on  the 
strength  of  the  shock  wave.  Thus  the  extent  of  amplification  should  decrease  in 
time  as  the  shock  wave  propagates  into  the  background  with  decreasing  intensity. 
Figure  34  a  shows  statistics  for  turbulent  kinetic  energy  at  t  —  1.9  2.25  fts 

and  2.5  ps  plotted  as  a  function  of  r  at  z  corresponding  to  the  center  of  the 
plasma.  Note  that  the  turbulence  levels  get  suppressed  in  the  plasma  core.  Also 
note  that  the  amplification  in  turbulence  levels  across  the  shock  wave  decreases 
in  time  due  to  decreasing  intensity  of  the  shock  wave.  Figure  34  b  shows  statistics 
for  divergence  in  time.  The  two  successive  peaks  indicate  expansion  in  the  core 
followed  by  compression  across  the  blast  wave.  Note  that  at  any  given  instant 
of  time  the  effect  of  compression  on  the  background  turbulence  is  much  more 
pronounced  compared  to  the  effect  of  expansion  in  the  plasma  core. 

4  Summary 

This  work  uses  direct  numerical  simulation  to  study  the  thermal  effect  of  a  re¬ 
gion  of  plasma  on  quiescent  air.  The  simulations  solve  the  compressible  ideal 
gas  Navier-Stokes  equations  using  Fourier  spectral  methods.  A  shock  capturing 
scheme  is  incorporated  to  account  for  the  strong  shock  waves.  Also,  a  logarith¬ 
mic  formulation  for  the  continuity  equation  is  developed  to  handle  low  densities 
at  the  core  of  the  plasma. 

Two  idealizations  of  the  plasma  have  been  considered  -  spherical  and  tear-drop 
shaped.  The  spherical  idealization  can  be  compared  to  classical  solution  for 
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Figure  35;  (a)  Radial  location  of  blast  wave  corresponding  to  Ms  =  1,05  for 
different  A T0)  (b)  Time  required  for  Ms  to  reduce  to  1.05  for  different  AT0. 


strong  shock  conditions.  The  post  energy  deposition  process  can  be  divided 
into  shock  formation  and  propagation  phases.  During  shock  formation,  large 
pressure  gradients  accelerate  the  flow  in  the  radial  direction.  This  results  in 
the  formation  of  a  shock  wave  whose  intensity  increases  in  time  until  it  reaches 
a  maximum.  The  shock  wave  then  propagates  radially*  its  strength  decreasing 
as  a  result.  Density  in  the  core  of  the  plasma  decreases  while  the  shock  forms 
and  increases  when  the  shock  propagates  radially.  This  behavior  of  density  is 
explained  by  conservation  of  mass.  The  shock  radius  and  velocity  at  the  shock 
front  are  compared  to  strong  shock  solutions  and  are  observed  to  compare  well 
during  initial  times  but  later  deviate  towards  the  acoustic  limit. 

The  tear-drop  shaped  idealization  resembles  the  initial  shape  of  the  plasma 
as  observed  in  experiments.  Shock  radius  and  jumps  in  fluid  properties  at  the 
shock  front  normal  to  the  axis  of  the  plasma  compare  reasonably  to  results  from 
experiments  and  other  simulations.  A  tear-drop  shaped  blast  wave  propagates 
down  the  domain,  becoming  increasingly  spherical  in  time.  The  pressure  gra¬ 
dients  are  observed  to  be  most  intense  normal  to  the  axis  of  the  plasma.  Also 
propagation  of  a  curved  shock  wave  through  the  domain  is  found  to  misalign 
the  density  and  pressure  gradients,  which  generates  baroclinic  vorticity  in  the 
flow.  Budgets  for  the  vorticity  transport  equation  show  that  vorticity  generated 
initially  through  baroclinic  means  eventually  decreases  through  bulk  expansion 
in  the  flow. 

For  control  applications,  it  is  useful  to  estimate  the  extent  to  which  the  energy 
addition  affects  the  flow.  The  blast  wave  that  is  generated  as  a  result  of  energy 
deposition,  decays  in  intensity  until  eventually  it  becomes  an  acoustic  wave.  In 
order  to  estimate  the  spatial  extent  to  which  the  blast  wave  affects  the  back¬ 
ground  flow,  we  (arbitrarily)  define  a  cut-off  value  of  A/5  =  LOS.  Beyond  this 
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M3  the  effect  of  the  blast  wave  is  considered  small.  Figure  35a  shows  the  radial 
location  of  the  blast  wave  when  A/a  decays  to  its  cut-off  value,  for  different  AT o- 
The  location  of  the  blast  wave  rc  is  normalized  with  respect  to  the  initial  size 
of  the  plasma  r0  =  n/ 20.  The  time  required  to  reach  this  location  tc  is  shown 
in  figure  35b.  tc  is  normalized  with  respect  to  f0  which  is  the  time  required  by 
an  acoustic  wave  in  the  background  medium  to  travel  a  distance  equal  to  the 
initial  size  of  the  plasma.  These  results  are  for  the  spherical  idealization  of  the 
plasma.  Note  that  when  ATq  varies  from  30  to  240,  rc/ro  varies  from  3.88  to 
9.8  while  tc/i0  varies  from  5.52  to  13.11. 
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5  Appendix  :  Eigen  vector  Matrix  formulation 
for  a  generalized  coordinate  system 


Let  the  Euler  equations  in  conservative  form  be  represented  by 

au 


at 


+  V  ■  F  =  0 


U  is  the  vector  of  the  conserved  variables 

f  9 

pu 

U  =  pv 
pw 

\  pe  0  / 


F  is  the  flux  vector  given  by 


F  = 


/  pVn  ^ 

puvn  +  pnx 
pvvn  4-  pny 
pwvn  +  pnz 

\  P^QVr,  / 


where 


vn  =  unx  +  vny  +  wnt 


And 


(50) 


(51) 


nx 2  +  ny2  +  nx2  =  1 
And  ho  is  the  total  enthalpy  given  by 

ho  —  h  +  Ck- 

Where  h  is  enthalpy  and  is  given  by 

h  —  e  +  p/p* 

Let  the  equation  of  state  for  the  gas  be  given  by 

P  =  pR(T)T 

Then  the  elemental  change  in  pressure  dp  can  be  written  as 
dp  =  RTdp  +  p{R  +  T—'jdT 
Now  since  e  =  e(T)  only 


(52) 
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Where  the  internal  energy  e  can  be  written  in  terms  of  the  conservative  variables 
Using  the  above  relations  the  generalized  Jacobian  matrix  can  be  computed 
as 


J  = 
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Where 

eg  —  e  +  ejt 

is  the  total  energy  and  ek  is  the  kinetic  energy  given  by 


ek  =  -utUi 


And  the  factor  A  is  given  by 


,dR\dT 
dTJde 

The  generalized  eigen  values  of  the  above  matrix  is  obtained  as  : 


/  d/?  \  i 

A-{R  +  Td fide 


A  =  (u„-ci,  vn,  w„+ci,  v„,  vn  ) 
Where  the  modified  speed  of  sound  c\  is  given  by 

Cl  =  «1  +  v4)/?r)1/2 

Let  a  new  set  of  variables  be  defined  such  that 
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Also  define  ek*  such  that 


e2  =  e  -  ej 
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e*'=e  fc-ea  (64) 

Using  these  new  variables  a  possible  set  of  the  right  eigen  vector  matrices  are 
obtained  as 
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Then  the  corresponding  set  of  left  eigen  vector  matrices  can  be  obtained  as 

(L]  =  [/?)-'  (65) 

and  they  are  given  by 
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